Analyses of seasonal and latitudinal distribution of larval fish in a changing western boundary current
Charles Hinchliffe
Summary
This script contains code to reproduce the analyses and figures from Hinchliffe C et al 2020, “Larval fish distribution in a changing western boundary current” (Link).
There are 3 major analyses within the manuscript. The first two involve analysis of the seasonal and latitudinal variation in larval fish abundance and taxa richness using generalised additive mixed-models (GAMMs) and can be found in section 2) and 3) in this script. The third major analyses involves investigation of the community composition of larval fish across the study area, first in relation to Sea Surface Temperature (SST), followed by Latitude (which appears in the supplementary material of the manuscript), and can be found in section 4) of this script. The figures associated with each analysis are at the end of each section, with there figure number as they appear in the manuscript and figure caption.
1) Data manipulation and cleaning
This section of the script cleans data and creates variables required for the subsequent analyses.
load required packages for all analyses and figure building.
library(tidyverse)
library(mgcv)
library(vegan)
library(car)
library(lubridate)
library(pracma)
#install_github('gavinsimpson/gratia@v0.2-3')
library(gratia)
library(lctools)
library(ggplot2)
library(gridGraphics)
library(cowplot)
library(ggpubr)
library(grid)
library(ggplotify)
library(mvabund)
#install ecoCopula
#devtools::install_github("gordy2x/ecoCopula",force=TRUE)
library(ecoCopula)
library(quantreg)bring in data
nimo.data <- read.csv("data/all_data.csv", header=TRUE)
#head(nimo.data)Make bathymetry absolute values in km units
nimo.data$Bathy <- (abs(nimo.data$Bathy)/1000)
#summary(nimo.data$Bathy)Create date variables
nimo.data$DateCaught <- as.Date(nimo.data$Date, "%d/%m/%Y")
nimo.data$Month <- month(nimo.data$Date)
#summary(nimo.data$Month)
nimo.data$Year <- year(nimo.data$DateCaught)
nimo.data$Day <- day(nimo.data$DateCaught)
#summary(nimo.data$Day)
nimo.data$dayofyear <- yday(nimo.data$DateCaught)
#summary(nimo.data$dayofyear)
#hist(nimo.data$dayofyear)Create season variables
nimo.data$Month <- as.factor(nimo.data$Month) #month as factor
#make season variable
nimo.data$Season <- NA_character_
nimo.data$Season[nimo.data$Month %in% c(1,2,12)] <- "summer"
nimo.data$Season[nimo.data$Month %in% 3:5] <- "autumn"
nimo.data$Season[nimo.data$Month %in% 6:8] <- "winter"
nimo.data$Season[nimo.data$Month %in% 9:11] <- "spring"
nimo.data$Season <- as.factor(nimo.data$Season)
nimo.data$Month <- as.numeric(nimo.data$Month) Create time period variables
nimo.data$Decade[nimo.data$Year <= 1999] <- "1990s"
nimo.data$Decade[nimo.data$Year <= 2009 & nimo.data$Year > 1999] <- "2000s"
nimo.data$Decade[nimo.data$Year <= 2020 & nimo.data$Year > 2009] <- "2010s"
nimo.data$Period[nimo.data$Year <= 1998] <- "pre1998"
nimo.data$Period[nimo.data$Year >= 1999] <- "post1998"
nimo.data$Decade <- as.factor(nimo.data$Decade)
#summary(nimo.data$Decade)
nimo.data$Period <- as.factor(nimo.data$Period)
#summary(nimo.data$Period)Create latitudinal region variable
nimo.data$Region[nimo.data$Latitude <= -40] <- "South"
nimo.data$Region[nimo.data$Latitude <= -35 & nimo.data$Latitude > -40] <- "Mid-South"
nimo.data$Region[nimo.data$Latitude <= -30 & nimo.data$Latitude > -35] <- "Mid-North"
nimo.data$Region[nimo.data$Latitude > -30] <- "North"
nimo.data$Region <- as.factor(nimo.data$Region)
nimo.data$Region <- as.factor(nimo.data$Region)
nimo.data$Region <- factor(nimo.data$Region, ordered = TRUE, levels = c("North", "Mid-North", "Mid-South","South"))
#summary(nimo.data$Region) Create mid-point standardised sampling depth variable
nimo.data$stand_depth <- as.character(nimo.data$Gear_depth_m)
nimo.data$stand_depth[nimo.data$stand_depth == "0-25"] <- "25" #double mid point
nimo.data$stand_depth[nimo.data$stand_depth == "0-30"] <- "30"
nimo.data$stand_depth[nimo.data$stand_depth == "0-35"] <- "35"
nimo.data$stand_depth[nimo.data$stand_depth == "0-40"] <- "40"
nimo.data$stand_depth[nimo.data$stand_depth == "0-50"] <- "50"
nimo.data$stand_depth[nimo.data$stand_depth == "0-80"] <- "80"
nimo.data$stand_depth[nimo.data$stand_depth == "25-50"] <- "75"
nimo.data$stand_depth[nimo.data$stand_depth == "50-75"] <- "125"
nimo.data$stand_depth[nimo.data$stand_depth == "75-100"] <- "175"
nimo.data$stand_depth <- as.factor(nimo.data$stand_depth)
nimo.data$stand_depth <- as.numeric(nimo.data$stand_depth)
nimo.data$stand_depth <- (nimo.data$stand_depth/2)Create sampling depth factors
nimo.data$depth_sect[nimo.data$stand_depth >= 0 & nimo.data$stand_depth < 11] <- "a) 0-10"
nimo.data$depth_sect[nimo.data$stand_depth >= 11 & nimo.data$stand_depth < 50] <- "b) 11-50"
nimo.data$depth_sect[nimo.data$stand_depth >= 51 & nimo.data$stand_depth < 100] <- "c) 51-100"
nimo.data$depth_sect[nimo.data$stand_depth >= 101 & nimo.data$stand_depth < 150] <- "d) 101-150"
nimo.data$depth_sect[nimo.data$stand_depth >= 151] <- "e) 151+"
nimo.data$depth_sect <- as.factor(nimo.data$depth_sect)Creating diversity variables
larval.spp<-nimo.data[1:3193,22:241]
abundance <- rowSums(larval.spp)
richness <- specnumber(larval.spp) # richness
shannon <- diversity(larval.spp) #shannon weaver index from vegan package
J <- shannon/log(specnumber(larval.spp)) # Evenness
shannon_effective <- exp(diversity(larval.spp))
nimo.data$shannon <- shannon
nimo.data$evenness <- J
nimo.data$richness <- richness
nimo.data$abundance <- abundance
nimo.data$shannon_effective <- shannon_effectiveSeparate Australian west coast and east coast
Samples into separate data frames
eastcoast <- subset(nimo.data, nimo.data$Longitude >= 130)
westcoast <- subset(nimo.data, nimo.data$Longitude <= 130)Check SST and Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST) correlation on east coast
eastcoast$SST <- as.numeric(eastcoast$SST)
eastcoast$HadISST <- as.numeric(eastcoast$HadISST)
cor.test(eastcoast$SST,eastcoast$HadISST, method = c("pearson", "kendall", "spearman"))##
## Pearson's product-moment correlation
##
## data: eastcoast$SST and eastcoast$HadISST
## t = 108.89, df = 1443, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9383103 0.9495273
## sample estimates:
## cor
## 0.9441919#plot(x = eastcoast$SST,y = eastcoast$HadISST)Separate inshore and offshore
Sites on the east coast based on position on continental shelf (boundary set at 200m bathymetry).
inshore <- subset(eastcoast, eastcoast$Bathy <= .200)
offshore <- subset(eastcoast, eastcoast$Bathy >= .200)2) Abundance GAMM
The following section contains code for the analysis of larval fish abundance using generalised additive mixed models (GAMM).
Testing linear correlation between predictors using Pearson correlation
preds <- data.frame(eastcoast$Latitude, eastcoast$HadISST, eastcoast$stand_depth, eastcoast$Bathy, eastcoast$dists_km, eastcoast$Temperature_C, eastcoast$SST)
cor(preds)## eastcoast.Latitude eastcoast.HadISST eastcoast.stand_depth eastcoast.Bathy eastcoast.dists_km eastcoast.Temperature_C eastcoast.SST
## eastcoast.Latitude 1.0000000 0.7939551 -0.20202890 NA 0.29394633 NA NA
## eastcoast.HadISST 0.7939551 1.0000000 -0.16558965 NA 0.27806479 NA NA
## eastcoast.stand_depth -0.2020289 -0.1655896 1.00000000 NA 0.00426594 NA NA
## eastcoast.Bathy NA NA NA 1 NA NA NA
## eastcoast.dists_km 0.2939463 0.2780648 0.00426594 NA 1.00000000 NA NA
## eastcoast.Temperature_C NA NA NA NA NA 1 NA
## eastcoast.SST NA NA NA NA NA NA 1Model selection
# step 1) Full model
model1 <- gam(abundance ~ Season + Period + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(dists_km, bs = "cs", k = 5) + s(HadISST, bs = "cs", k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE, Select=TRUE)
concurvity(model1, full=T) #checking concurvity## para s(stand_depth) s(dists_km) s(HadISST) s(Latitude):Seasonautumn s(Latitude):Seasonspring s(Latitude):Seasonsummer s(Latitude):Seasonwinter s(Project_name)
## worst 1 0.3174216 0.8765001 0.9479951 0.9191817 0.9259444 0.9537419 0.9044561 1.0000000
## observed 1 0.1266766 0.7673233 0.9215439 0.8004886 0.8111269 0.9347195 0.7972925 0.9589919
## estimate 1 0.1626038 0.6044313 0.8125066 0.7369388 0.7975847 0.8583453 0.8191350 0.8065563# step 2) s(dists_km, bs = "cs", k = 5) + s(HadISST, bs = "cs", k = 5) removed due to concurvity
model2 <- gam(abundance ~ Season + Period + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE, Select=TRUE)
# Step 3) Period removed based on AIC score, checked with chi-sq test
model2i <- gam(abundance ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE, Select=TRUE)
AIC(model2,model2i) ## AIC reduced with period removed## df AIC
## model2 30.69146 32233.64
## model2i 29.84881 32231.60anova(model2,model2i, test = "Chisq") #non sig therefore can remove period (adds nothing to model)## Analysis of Deviance Table
##
## Model 1: abundance ~ Season + Period + exp(-Bathy) + s(stand_depth, bs = "cs",
## k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season,
## id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re")
## Model 2: abundance ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs",
## k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season,
## id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re")
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 2968.1 3580.7
## 2 2969.0 3580.7 -0.93486 0.021107# Step 4) linear Season term removed, but increases AIC and is sig dif in chi-sq test there retain Season in model
model3 <- gam(abundance ~ exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE, Select=TRUE)
AIC(model3,model2i) ## AIC increased when season removed. ## model2i best model here. ## df AIC
## model3 29.40232 32256.84
## model2i 29.84881 32231.60anova(model3,model2i, test = "Chisq") #sig dif therefore dont remove season## Analysis of Deviance Table
##
## Model 1: abundance ~ exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) +
## s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) +
## offset(log(Volume_m3)) + s(Project_name, bs = "re")
## Model 2: abundance ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs",
## k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season,
## id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re")
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 2969.4 3582.7
## 2 2969.0 3580.7 0.41619 1.9727 0.05599 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## therefore best model for abundance appears to be model2i, removing Period based on AIC and s(dists_km, bs = "cs", k = 5) + s(HadISST, bs = "cs", k = 5) based on concurvity.
#Run with method = REML for estimates
model2i <- gam(abundance ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "REML", fit=TRUE, Select=TRUE)
#Best model summary
summary(model2i) ##
## Family: Negative Binomial(0.586)
## Link function: log
##
## Formula:
## abundance ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs",
## k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season,
## id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.93971 0.30561 -6.347 2.19e-10 ***
## Seasonspring -0.55675 0.11191 -4.975 6.52e-07 ***
## Seasonsummer 0.11398 0.06852 1.664 0.09621 .
## Seasonwinter -0.31556 0.10983 -2.873 0.00406 **
## exp(-Bathy) 0.96264 0.13170 7.309 2.69e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(stand_depth) 3.866 4 582.6 < 2e-16 ***
## s(Latitude):Seasonautumn 2.377 4 303.0 < 2e-16 ***
## s(Latitude):Seasonspring 1.967 4 442.4 3.92e-08 ***
## s(Latitude):Seasonsummer 2.469 4 817.2 5.58e-10 ***
## s(Latitude):Seasonwinter 2.064 4 273.9 1.96e-08 ***
## s(Project_name) 9.487 10 230.2 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = -0.079 Deviance explained = 30.2%
## -REML = 16152 Scale est. = 1 n = 3000#par(mfcol=c(3,2))
#plot(model2i, shade = TRUE)
#checking concurvity
#concurvity(model2i,full=TRUE)
#checking assumptions
appraise(model2i)
#over disperison test
E1 <- resid(model2i, type = "pearson") #ok - less than 2
N <- nrow(eastcoast)
p <- length(coef(model2i))
Overdispersion <- sum(E1^2) / (N - p)
#paste('Overdispersion: ',Overdispersion) #"Overdispersion: 1.58845886919077"Checking spatial auto correlation of Abundance GAMM with Moran’s I If Moran’s I is ~0 then there is no spatial autocorrelation.
res <- residuals(model2i)
#res
Coords <-cbind(eastcoast$Latitude, eastcoast$Longitude)
bw <- 6
mI <-moransI(Coords,bw,res)
moran.table <-matrix(data=NA,nrow=1,ncol=6)
col.names <-c("Moran's I", "Expected I", "Z resampling", "P-value resampling","Z randomization", "P-value randomization")
colnames(moran.table) <- col.names
moran.table[1,1] <- mI$Morans.I
moran.table[1,2] <- mI$Expected.I
moran.table[1,3] <- mI$z.resampling
moran.table[1,4] <- mI$p.value.resampling
moran.table[1,5] <- mI$z.randomization
moran.table[1,6] <- mI$p.value.randomization
print(moran.table)Plotting output
Plotting effects of the larval fish abundance GAMM.
Figure 2.
Figure 2. Estimated terms in GAMM model of total larval abundance, where the black line is equal to the mean effect and the grey ribbon is the 95% Confidence Interval. A) Partial linear effect of exp(-Bathymetry) in kilometres. B-D) Estimated smoothers for sampling depth (meters), latitude during Summer, latitude during Autumn, latitude during Winter, and latitude during Spring. n = 3006.
# plot Bathymetry effect
z <- evaluate_parametric_term(model2i, "exp(-Bathy)")
#z
pbath <- ggplot(z, aes(x = value, y = partial)) +
geom_line() +
geom_ribbon(data=z,aes(ymin=upper,ymax=lower),alpha=0.3) +
xlab("exp(-Bathy)") +
ylab("Partial effect of exp(-Bathy)") +
geom_rug(sides ="b", color="black") +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(plot.margin = margin(30,20,32,24))
#plot smooth effects and save as grob file
p1 <- as.grob(function() plot(model2i, shade = TRUE, select = 1))
#p1
p2 <- as.grob(function() plot(model2i, shade = TRUE, select = 2))
#p2
p3 <- as.grob(function() plot(model2i, shade = TRUE, select = 3))
#p3
p4 <- as.grob(function() plot(model2i, shade = TRUE, select = 4))
#p4
p5 <- as.grob(function() plot(model2i, shade = TRUE, select = 5))
#p5
figure_2 <- plot_grid(pbath, p1, p4, p2, p5, p3, labels = c('A', 'B', 'C', 'D', 'E', 'F'), label_size = 12, nrow = 3, ncol = 2)
figure_2
3) Richness GAMM
The following section contains code for the analysis of larval fish taxa richness using generalised additive mixed models (GAMM).
Model selection
# step 1) Full model
model1 <- gam(richness ~ Season + Period + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(dists_km, bs = "cs", k = 5) + s(HadISST, bs = "cs", k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE, Select=TRUE)
#concurvity(model1, full=T) #checking concurvity
# step 2) s(dists_km, bs = "cs", k = 5) + s(HadISST, bs = "cs", k = 5) removed due to concurvity
model2 <- gam(richness ~ Season + Period + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE, Select=TRUE)
# Step 3) Period removed based on AIC score, checked with chi-sq test
model2i <- gam(richness ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE, Select=TRUE)
AIC(model2,model2i) ## AIC reduced with period removed## df AIC
## model2 30.49724 19565.42
## model2i 29.82409 19561.91#anova(model2,model2i, test = "Chisq") #non sig therefore can remove period (adds nothing to model)
# Step 4) linear Season term removed, but increases AIC and is sig dif in chi-sq test there retain Season in model
model3 <- gam(richness ~ exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE, Select=TRUE)
AIC(model3,model2i) ## AIC increased when season removed. ## model2i best model here. ## df AIC
## model3 30.68741 19623.27
## model2i 29.82409 19561.91anova(model3,model2i, test = "Chisq") #sig dif therefore don't remove season## Analysis of Deviance Table
##
## Model 1: richness ~ exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(Latitude,
## bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) +
## s(Project_name, bs = "re")
## Model 2: richness ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) +
## s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) +
## offset(log(Volume_m3)) + s(Project_name, bs = "re")
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 2968.4 3360.4
## 2 2969.1 3358.8 -0.66848 1.5397## therefore best model for richness appears to be model2i, removing Period based on AIC and s(dists_km, bs = "cs", k = 5) + s(HadISST, bs = "cs", k = 5) based on concurvity.
# Run with method = REML for estimates
model2i <- gam(richness ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "REML", fit=TRUE, Select=TRUE)
#Best model summary
summary(model2i) ##
## Family: Negative Binomial(2.032)
## Link function: log
##
## Formula:
## richness ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) +
## s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) +
## offset(log(Volume_m3)) + s(Project_name, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.98100 0.20233 -19.676 < 2e-16 ***
## Seasonspring -0.45834 0.07145 -6.415 1.41e-10 ***
## Seasonsummer 0.11461 0.04220 2.716 0.00661 **
## Seasonwinter -0.37823 0.06931 -5.457 4.85e-08 ***
## exp(-Bathy) 0.83862 0.07820 10.724 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(stand_depth) 3.866 4 604.4 < 2e-16 ***
## s(Latitude):Seasonautumn 2.573 4 483.0 < 2e-16 ***
## s(Latitude):Seasonspring 2.169 4 749.8 7.61e-09 ***
## s(Latitude):Seasonsummer 2.714 4 178.0 0.0128 *
## s(Latitude):Seasonwinter 2.233 4 435.8 1.76e-13 ***
## s(Project_name) 9.557 10 570.9 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = -0.0346 Deviance explained = 45.4%
## -REML = 9823.8 Scale est. = 1 n = 3000#par(mfcol=c(3,2))
#plot(model2i, shade = TRUE) ### spline plots
#checking concurvity
#concurvity(model2i,full=TRUE)
#checking assumptions
appraise(model2i)
#overdisperison test
E1 <- resid(model2i, type = "pearson") #ok - less than 2
N <- nrow(eastcoast)
p <- length(coef(model2i))
Overdispersion <- sum(E1^2) / (N - p)
#paste('Overdispersion: ',Overdispersion) #"Overdispersion: 1.01888080182827"Checking spatial auto correlation of Taxa Richness GAMM with Moran’s I If Moran’s I is ~0 then there is no spatial autocorrelation.
res <- residuals(model2i)
#res
Coords <-cbind(eastcoast$Latitude, eastcoast$Longitude)
bw <- 6
mI <-moransI(Coords,bw,res)
moran.table <-matrix(data=NA,nrow=1,ncol=6)
col.names <-c("Moran's I", "Expected I", "Z resampling", "P-value resampling","Z randomization", "P-value randomization")
colnames(moran.table) <- col.names
moran.table[1,1] <- mI$Morans.I
moran.table[1,2] <- mI$Expected.I
moran.table[1,3] <- mI$z.resampling
moran.table[1,4] <- mI$p.value.resampling
moran.table[1,5] <- mI$z.randomization
moran.table[1,6] <- mI$p.value.randomization
print(moran.table)Plotting output
Plotting effects of the larval fish taxa richness GAMM.
Figure 3.
Figure 3. Estimated terms in GAMM model of larval taxa richness, where the black line is equal to the mean effect and the grey ribbon is 95% Confidence Interval. A) Partial linear effect of exp(-Bathymetry) in kilometres. B) Non-linear effect of sampling depth (meters). C) Non-linear effect of latitude during Summer. D) Non-linear effect of latitude during Autumn. E) Non-linear effect of latitude during Winter. F) Non-linear effect of latitude during Spring. n = 3006.
z <- evaluate_parametric_term(model2i, "exp(-Bathy)")
#z
pbath <- ggplot(z, aes(x = value, y = partial)) +
geom_line() +
geom_ribbon(data=z,aes(ymin=upper,ymax=lower),alpha=0.3) +
xlab("exp(-Bathy)") +
ylab("Partial effect of exp(-Bathy)") +
geom_rug(sides ="b", color="black") +
theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(plot.margin = margin(30,20,32,24))
p1 <- as.grob(function() plot(model2i, shade = TRUE, select = 1))
#p1
p2 <- as.grob(function() plot(model2i, shade = TRUE, select = 2))
#p2
p3 <- as.grob(function() plot(model2i, shade = TRUE, select = 3))
#p3
p4 <- as.grob(function() plot(model2i, shade = TRUE, select = 4))
#p4
p5 <- as.grob(function() plot(model2i, shade = TRUE, select = 5))
#p5
figure_3 <- plot_grid(pbath, p1, p4, p2, p5, p3, labels = c('A', 'B', 'C', 'D', 'E', 'F'), label_size = 12, nrow = 3, ncol = 2)
figure_3
4) Multivariate Analyses
This section of the scripts uses an intercept-only multivariate generalised linear model (MGLM) model of the larval fish community in a Gaussian copula graphical model (GCGM), and quantile regression of the , to observe the effects of temperature (Figure 4), latitude (Supplementary Figure 2) and time period (Supplementary Figure 3).
Select inshore species in matrix to analyse
elements <- inshore[1:2428,22:239]
myvars <- names(elements) %in% c("Acropomatidae_Synagrops.spp_37311949",
"Aploactinidae_other_37290000",
"Bothidae_Crossorhombus.spp_37460907",
"Bovichtidae_Pseudaphritis.urvilli_37403003",
"Callanthiidae_other_37311957",
"Cepolidae_Acanthocepola.spp_37380901",
"Cepolidae_Owstonia.spp_37380903",
"Cepolidae_other_37380000",
"Cetomimidae_37132000",
"Chandidae_other_37310900",
"Engraulidae_other_37086000",
"Ipnopidae_37123000",
"Latridae_other_37378000",
"Leptobramidae_37357905",
"Microcanthidae_other_37361900",
"Ophidiidae_Brotula.spp_37228912",
"Plesiopidae_37316000",
"Terapontidae_other_37321000",
"Xiphiidae_Xiphias.gladius_37442001")
elements <- elements[!myvars] #removing species which never occur on in less than 200m bathymetry on the east coast.
elements <- mvabund(elements[,])Fitting GCGM with intercept only MGLM
#boxplot(elements) # inspect data ranges
#meanvar.plot(elements) # check mean variance relationship
#fit <- manyglm(elements ~Latitude, family = "negative.binomial", offset = log(Volume_m3),data=inshore)
fit_null <- manyglm(elements ~1, family = "negative.binomial", offset = log(Volume_m3),data=inshore)
fit_cop <- cord(fit_null)
plot(fit_cop,biplot=TRUE)
score_1 <- fit_cop$scores[[1]][1,]
score_2 <- fit_cop$scores[[1]][2,]Plot GCGM bivariate scores 1 and 2 with SST, and then latitude.
plot_dat <- data.frame(score_1,score_2,inshore$Region, inshore$Latitude, inshore$HadISST, inshore$Period)
SST_bivariate_plot <- ggplot(plot_dat,aes(x=score_1,y=score_2,col=inshore.HadISST))+
geom_point() +
xlab("Score 1") +
ylab("Score 2") +
scale_color_continuous(high = "red", low = "blue") +
theme_bw() +
theme(legend.position=c(0.85, .9), legend.title = element_text(size=10, face="bold"),
legend.text = element_text(size = 10, face = "bold"))
LAT_bivariate_plot <- ggplot(plot_dat,aes(x=score_1,y=score_2,col=inshore.Latitude))+
geom_point() +
xlab("Score 1") +
ylab("Score 2") +
scale_color_continuous(high = "red", low = "blue") +
theme_bw() +
theme(legend.position=c(0.85, .9), legend.title = element_text(size=10, face="bold"),
legend.text = element_text(size = 10, face = "bold"))Quantile regression of bivariate scores against SST
Quantile regression of score 1 from GCGM against SST
qs <- c(0.05, 0.5, 0.95) ### What quantiles are we interested in... 5%, 50% and 95%
#quantile regression for score 1
qr1 <- rq(score_1 ~ inshore.HadISST, data=plot_dat, tau = qs)
SM <- summary(qr1, se = "boot", bsmethod = "wild", R = 1000) ## method of calculating SE and bootstrap method if applicable
SM##
## Call: rq(formula = score_1 ~ inshore.HadISST, tau = qs, data = plot_dat)
##
## tau: [1] 0.05
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -1.46418 0.08469 -17.28774 0.00000
## inshore.HadISST 0.02291 0.00444 5.15674 0.00000
##
## Call: rq(formula = score_1 ~ inshore.HadISST, tau = qs, data = plot_dat)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -1.50302 0.09426 -15.94554 0.00000
## inshore.HadISST 0.06436 0.00526 12.23176 0.00000
##
## Call: rq(formula = score_1 ~ inshore.HadISST, tau = qs, data = plot_dat)
##
## tau: [1] 0.95
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -4.01042 0.53902 -7.44027 0.00000
## inshore.HadISST 0.29237 0.02887 10.12823 0.00000Plotting quantile regression score 1 from GCGM against SST
### Build a ggplot
p <- ggplot(plot_dat, aes(x = inshore.HadISST, y = score_1)) +
geom_point(size = 3, na.rm = TRUE, shape = 1)
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + geom_abline(intercept = qr1$coefficients[1,1], slope = qr1$coefficients[2,1],
linetype = 'longdash', col = "blue") # 5th quantile
p <- p + geom_abline(intercept = qr1$coefficients[1,2], slope = qr1$coefficients[2,2], col = "blue")
# linetype = 'longdash' ) # 50th quantile
p <- p + geom_abline(intercept = qr1$coefficients[1,3], slope = qr1$coefficients[2,3],
linetype = 'longdash', col = "blue") # 95th quantile
p <- p + xlab("SST (°C)")
p <- p + ylab("Score 1")
p <- p + theme(axis.text=element_text(size=16, face = "bold"))
p <- p + theme(axis.title=element_text(size=16, face = "bold"))
p_score1 <- p + theme(legend.position=c(0.15, .8), legend.title = element_text(size=16, face="bold"),
legend.text = element_text(size = 16, face = "bold"))
#p_score1Quantile regression of score 2 from GCGM against SST
qs <- c(0.05, 0.5, 0.95) ### What quantiles are we interested in... 5%, 50% and 95%
#quantile regression for score 1
qr2 <- rq(score_2 ~ inshore.HadISST, data=plot_dat, tau = qs)
SM <- summary(qr2, se = "boot", bsmethod = "wild", R = 1000) ## method of calculating SE and bootstrap method if applicable
SM##
## Call: rq(formula = score_2 ~ inshore.HadISST, tau = qs, data = plot_dat)
##
## tau: [1] 0.05
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -1.84492 0.13320 -13.85082 0.00000
## inshore.HadISST 0.04684 0.00723 6.47927 0.00000
##
## Call: rq(formula = score_2 ~ inshore.HadISST, tau = qs, data = plot_dat)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -2.17395 0.12361 -17.58653 0.00000
## inshore.HadISST 0.10387 0.00657 15.80505 0.00000
##
## Call: rq(formula = score_2 ~ inshore.HadISST, tau = qs, data = plot_dat)
##
## tau: [1] 0.95
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -4.28122 0.13235 -32.34750 0.00000
## inshore.HadISST 0.28888 0.00879 32.86453 0.00000Plotting quantile regression score 2 from GCGM against SST
### Build a ggplot
p <- ggplot(plot_dat, aes(x = inshore.HadISST, y = score_2)) +
geom_point(size = 3, na.rm = TRUE, shape = 1)
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) ## set up plot area
p <- p + geom_abline(intercept = qr2$coefficients[1,1], slope = qr2$coefficients[2,1],
linetype = 'longdash', col = "blue") # 5th quantile
p <- p + geom_abline(intercept = qr2$coefficients[1,2], slope = qr2$coefficients[2,2], col = "blue")
# linetype = 'longdash' ) # 50th quantile
p <- p + geom_abline(intercept = qr2$coefficients[1,3], slope = qr2$coefficients[2,3],
linetype = 'longdash', col = "blue") # 95th quantile
p <- p + xlab("SST (°C)")
p <- p + ylab("Score 2")
p <- p + theme(axis.text=element_text(size=16, face = "bold"))
p <- p + theme(axis.title=element_text(size=16, face = "bold"))
p_score2 <- p + theme(legend.position=c(0.15, .8), legend.title = element_text(size=16, face="bold"),
legend.text = element_text(size = 16, face = "bold"))
#p_score2Plotting the GCGM bivariate scores 1 and 2 with SST, and quantile regression plots, as per Figure 4 in the manuscript.
Figure 4.
Figure 4. A) GCGM bivariate scores 1 and 2 from the MGLM intercept-only model of all larval fish taxa abundances . Each point represents an assemblage and the colour of the points represents SST in °C. B) GCGM bivariate score 1 by SST and C) GCGM bivariate score 2 by SST. Blue lines represent 0.05, 0.5 and 0.95 quantiles. n = 2428.
figure_4 <- plot_grid(SST_bivariate_plot, p_score1, p_score2, labels = c('A', 'B', 'C'), label_size = 12, nrow = 1, ncol = 3)
figure_4
Univariate tests for MGLM with SST
The univariate tests for the MGLM with SST were run on the University of New South Wales computational cluster ‘Katana’ supported by Research Technology Services. The r script included in the repo ‘manyglm_uni_test.R’ contains the analyses. Here, we import the output to observe the results.
Read in rds and display output
Uni_test_manyglm_output <- readRDS("outputs/Uni_test_manyglm_output.rds")
#Uni_test_manyglm_output
#Uni_test_manyglm_output$uni.p
uni_p <- as.data.frame(Uni_test_manyglm_output$uni.p)
#uni_p
uni_ts <- as.data.frame(Uni_test_manyglm_output$uni.test)
#uni_ts
useful_output <- rbind(uni_p, uni_ts)
#useful_output
useful_output <- t(useful_output)
useful_output <- as.data.frame(useful_output)
sig_out <- subset(useful_output, useful_output$HadISST <= 0.05)
#sig_out
#sum(useful_output$HadISST1) #whole test dev
useful_output$dev_exp <- (useful_output$HadISST1/sum(useful_output$HadISST1))*100
#useful_output
useful_output <- useful_output %>%
dplyr::rename(
P = HadISST,
TS = HadISST1
) %>%
dplyr::select(P, TS, dev_exp)
head(useful_output)## P TS dev_exp
## Acanthuridae_37437000 0.0001 156.978730 1.1467679622
## Acropomatidae_Apogonops.anomalus_37311053 0.0001 95.979510 0.7011537639
## Acropomatidae_other_37311956 1.0000 0.111148 0.0008119635
## Ammodytidae_37425000 0.0002 56.514417 0.4128516210
## Anguilliformes_37990019 0.0001 271.926545 1.9864898281
## Antennariidae_37210000 0.4399 9.786945 0.0714960249#write.csv(useful_output, "manyglm_temp_uni.csv")Quantile regression of bivariate scores against Latitude
Now running quantile regression of score 1 from GCGM against Latitude
qs <- c(0.05, 0.5, 0.95) ### What quantiles are we interested in... 5%, 50% and 95%
#quantile regression for score 1
qr1 <- rq(score_1 ~ inshore.Latitude, data=plot_dat, tau = qs)
SM <- summary(qr1, se = "boot", bsmethod = "wild", R = 1000) ## method of calculating SE and bootstrap method if applicable
SM##
## Call: rq(formula = score_1 ~ inshore.Latitude, tau = qs, data = plot_dat)
##
## tau: [1] 0.05
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -0.46350 0.17303 -2.67868 0.00744
## inshore.Latitude 0.01580 0.00468 3.37658 0.00075
##
## Call: rq(formula = score_1 ~ inshore.Latitude, tau = qs, data = plot_dat)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 1.42931 0.11535 12.39081 0.00000
## inshore.Latitude 0.04724 0.00306 15.45606 0.00000
##
## Call: rq(formula = score_1 ~ inshore.Latitude, tau = qs, data = plot_dat)
##
## tau: [1] 0.95
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 9.58603 0.40829 23.47826 0.00000
## inshore.Latitude 0.22208 0.01057 21.01611 0.00000Plotting quantile regression score 1 from GCGM against Latitude
### Build a ggplot
p <- ggplot(plot_dat, aes(x = inshore.Latitude, y = score_1)) +
geom_point(size = 3, na.rm = TRUE, shape = 1)
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + geom_abline(intercept = qr1$coefficients[1,1], slope = qr1$coefficients[2,1],
linetype = 'longdash', col = "blue") # 5th quantile
p <- p + geom_abline(intercept = qr1$coefficients[1,2], slope = qr1$coefficients[2,2], col = "blue")
# linetype = 'longdash' ) # 50th quantile
p <- p + geom_abline(intercept = qr1$coefficients[1,3], slope = qr1$coefficients[2,3],
linetype = 'longdash', col = "blue") # 95th quantile
p <- p + xlab("Latitude (°)")
p <- p + ylab("Score 1")
p <- p + theme(axis.text=element_text(size=16, face = "bold"))
p <- p + theme(axis.title=element_text(size=16, face = "bold"))
p_score1 <- p + theme(legend.position=c(0.15, .8), legend.title = element_text(size=16, face="bold"),
legend.text = element_text(size = 16, face = "bold"))
#p_score1Now running quantile regression of score 1 from GCGM against Latitude
qs <- c(0.05, 0.5, 0.95) ### What quantiles are we interested in... 5%, 50% and 95%
#quantile regression for score 1
qr2 <- rq(score_2 ~ inshore.Latitude, data=plot_dat, tau = qs)
SM <- summary(qr2, se = "boot", bsmethod = "wild", R = 1000) ## method of calculating SE and bootstrap method if applicable
SM##
## Call: rq(formula = score_2 ~ inshore.Latitude, tau = qs, data = plot_dat)
##
## tau: [1] 0.05
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -0.08068 0.21772 -0.37058 0.71098
## inshore.Latitude 0.02462 0.00601 4.09674 0.00004
##
## Call: rq(formula = score_2 ~ inshore.Latitude, tau = qs, data = plot_dat)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 1.58956 0.16467 9.65318 0.00000
## inshore.Latitude 0.04990 0.00438 11.38512 0.00000
##
## Call: rq(formula = score_2 ~ inshore.Latitude, tau = qs, data = plot_dat)
##
## tau: [1] 0.95
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 8.61369 0.25249 34.11502 0.00000
## inshore.Latitude 0.19815 0.00614 32.26831 0.00000Plotting quantile regression score 2 from GCGM against Latitude
### Buiild a ggplot
p <- ggplot(plot_dat, aes(x = inshore.Latitude, y = score_2)) +
geom_point(size = 3, na.rm = TRUE, shape = 1)
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) ## set up plot area
p <- p + geom_abline(intercept = qr2$coefficients[1,1], slope = qr2$coefficients[2,1],
linetype = 'longdash', col = "blue") # 5th quantile
p <- p + geom_abline(intercept = qr2$coefficients[1,2], slope = qr2$coefficients[2,2], col = "blue")
# linetype = 'longdash' ) # 50th quantile
p <- p + geom_abline(intercept = qr2$coefficients[1,3], slope = qr2$coefficients[2,3],
linetype = 'longdash', col = "blue") # 95th quantile
p <- p + xlab("Latitude (°)")
p <- p + ylab("Score 2")
p <- p + theme(axis.text=element_text(size=16, face = "bold"))
p <- p + theme(axis.title=element_text(size=16, face = "bold"))
p_score2 <- p + theme(legend.position=c(0.15, .8), legend.title = element_text(size=16, face="bold"),
legend.text = element_text(size = 16, face = "bold"))
#p_score2Plotting the GCGM bivariate scores 1 and 2 with Latitude, and quantile regression plots, as per Supplementary Figure 2 in the manuscript.
Supp Figure 2.
Supplementary Figure 2. A) GCGM bivariate scores 1 and 2. Each point represents an assemblage and the colour of the points represents latitude. B) GCGM bivariate score 1 by latitude and C) GCGM bivariate score 2 by latitude. Blue lines represent 0.05, 0.5 and 0.95 quantiles. n = 2428.
supp_figure_2 <- plot_grid(LAT_bivariate_plot, p_score1, p_score2, labels = c('A', 'B', 'C'), label_size = 12, nrow = 1, ncol = 3)
supp_figure_2
Observing difference in community composition of larval fish taxa with latitude before and after 1998 (coinciding with the 1998/99 El Nino event)in the bivariate scores of the GCGM with intercept only MGLM.
Figure 5.
Figure 5. A) GCGM bivariate scores 1 and 2. Each point represents an assemblage and the colour of the points represents the ‘Period’, black circles are pre-1998 and orange circles are post-1998. B) GCGM bivariate score 1 by latitude and C) GCGM bivariate score 2 by latitude. Trend lines are fit with a loess smoother.
my_cols <- c("#E69F00", "#000000")
p1<- ggplot(plot_dat,aes(x=score_1,y=score_2))+
geom_point(aes(color = inshore.Period), size = 4, alpha = 0.5) +
scale_color_manual(values = my_cols) +
theme_classic() +
theme(legend.position=c(0.85, .9), legend.title = element_text(size=10, face="bold"),
legend.text = element_text(size = 10, face = "bold"))
#p1
p2 <- ggplot(plot_dat,aes(x=inshore.Latitude, y=score_1, col=inshore.Period))+
geom_point(aes(color = inshore.Period), size = 4, alpha = 0.5, shape = 1) +
scale_color_manual(values = my_cols) +
theme_classic() +
xlab("Latitude (°C)") +
theme(legend.position="none") +
geom_smooth(method = "loess")
p3 <- ggplot(plot_dat,aes(x=inshore.Latitude, y=score_2, col=inshore.Period))+
geom_point(aes(color = inshore.Period), size = 4, alpha = 0.5, shape = 1) +
scale_color_manual(values = my_cols) +
theme_classic() +
xlab("Latitude (°C)") +
theme(legend.position="none") +
geom_smooth(method = "loess")
supp_figure_3 <- plot_grid(p1, p2, p3, labels = c('A', 'B', 'C'), label_size = 12, nrow = 1, ncol = 3)
supp_figure_3
Ssssion Info
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## setting value
## version R version 3.6.1 (2019-07-05)
## os macOS Mojave 10.14.6
## system x86_64, darwin15.6.0
## ui X11
## language (EN)
## collate en_AU.UTF-8
## ctype en_AU.UTF-8
## tz Australia/Sydney
## date 2020-07-07
##
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 3.6.0)
## acepack 1.4.1 2016-10-29 [1] CRAN (R 3.6.0)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0)
## backports 1.1.5 2019-10-02 [1] CRAN (R 3.6.0)
## base64enc 0.1-3 2015-07-28 [1] CRAN (R 3.6.0)
## BiocManager 1.30.9 2019-10-23 [1] CRAN (R 3.6.1)
## broom 0.5.2 2019-04-07 [1] CRAN (R 3.6.0)
## callr 3.3.2 2019-09-22 [1] CRAN (R 3.6.0)
## car * 3.0-5 2019-11-15 [1] CRAN (R 3.6.0)
## carData * 3.0-3 2019-11-16 [1] CRAN (R 3.6.0)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 3.6.0)
## checkmate 1.9.4 2019-07-04 [1] CRAN (R 3.6.0)
## cli 2.0.2 2020-02-28 [1] CRAN (R 3.6.0)
## cluster 2.1.0 2019-06-19 [1] CRAN (R 3.6.1)
## colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.6.0)
## cowplot * 1.0.0 2019-07-11 [1] CRAN (R 3.6.0)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.0)
## curl 4.3 2019-12-02 [1] CRAN (R 3.6.0)
## data.table 1.12.8 2019-12-09 [1] CRAN (R 3.6.0)
## DBI 1.1.0 2019-12-15 [1] CRAN (R 3.6.0)
## dbplyr 1.4.2 2019-06-17 [1] CRAN (R 3.6.0)
## desc 1.2.0 2018-05-01 [1] CRAN (R 3.6.0)
## devtools 2.2.1 2019-09-24 [1] CRAN (R 3.6.0)
## digest 0.6.25 2020-02-23 [1] CRAN (R 3.6.0)
## dplyr * 1.0.0 2020-05-29 [1] CRAN (R 3.6.2)
## ecoCopula * 0.0.0.9000 2020-07-03 [1] Github (gordy2x/ecoCopula@ea048e1)
## ellipsis 0.3.0 2019-09-20 [1] CRAN (R 3.6.0)
## evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.0)
## fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.0)
## farver 2.0.1 2019-11-13 [1] CRAN (R 3.6.0)
## forcats * 0.4.0 2019-02-17 [1] CRAN (R 3.6.0)
## foreign 0.8-72 2019-08-02 [1] CRAN (R 3.6.0)
## Formula 1.2-3 2018-05-03 [1] CRAN (R 3.6.0)
## fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.0)
## gdata 2.18.0 2017-06-06 [1] CRAN (R 3.6.0)
## generics 0.0.2 2018-11-29 [1] CRAN (R 3.6.0)
## ggplot2 * 3.2.1 2019-08-10 [1] CRAN (R 3.6.0)
## ggplotify * 0.0.5 2020-03-12 [1] CRAN (R 3.6.0)
## ggpubr * 0.2.3 2019-09-03 [1] CRAN (R 3.6.0)
## ggsignif 0.6.0 2019-08-08 [1] CRAN (R 3.6.0)
## glue 1.4.0 2020-04-03 [1] CRAN (R 3.6.1)
## gratia * 0.2-3 2020-07-01 [1] Github (gavinsimpson/gratia@245ac45)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 3.6.0)
## gridGraphics * 0.5-0 2020-02-25 [1] CRAN (R 3.6.0)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.0)
## gtools 3.8.1 2018-06-26 [1] CRAN (R 3.6.0)
## haven 2.2.0 2019-11-08 [1] CRAN (R 3.6.0)
## Hmisc 4.4-0 2020-03-23 [1] CRAN (R 3.6.0)
## hms 0.5.3 2020-01-08 [1] CRAN (R 3.6.0)
## htmlTable 1.13.3 2019-12-04 [1] CRAN (R 3.6.0)
## htmltools 0.4.0 2019-10-04 [1] CRAN (R 3.6.0)
## htmlwidgets 1.5.1 2019-10-08 [1] CRAN (R 3.6.0)
## httr 1.4.1 2019-08-05 [1] CRAN (R 3.6.0)
## jpeg 0.1-8 2014-01-23 [1] CRAN (R 3.6.0)
## jsonlite 1.6.1 2020-02-02 [1] CRAN (R 3.6.0)
## knitr 1.26 2019-11-12 [1] CRAN (R 3.6.0)
## labeling 0.3 2014-08-23 [1] CRAN (R 3.6.0)
## lattice * 0.20-38 2018-11-04 [1] CRAN (R 3.6.1)
## latticeExtra 0.6-29 2019-12-19 [1] CRAN (R 3.6.0)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.6.0)
## lctools * 0.2-8 2020-04-06 [1] CRAN (R 3.6.2)
## lifecycle 0.2.0 2020-03-06 [1] CRAN (R 3.6.0)
## lubridate * 1.7.4 2018-04-11 [1] CRAN (R 3.6.0)
## magrittr * 1.5 2014-11-22 [1] CRAN (R 3.6.0)
## MASS 7.3-51.4 2019-03-31 [1] CRAN (R 3.6.1)
## Matrix 1.2-17 2019-03-22 [1] CRAN (R 3.6.1)
## MatrixModels 0.4-1 2015-08-22 [1] CRAN (R 3.6.0)
## memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.0)
## mgcv * 1.8-28 2019-03-21 [1] CRAN (R 3.6.1)
## mice 3.8.0 2020-02-21 [1] CRAN (R 3.6.0)
## modelr 0.1.5 2019-08-08 [1] CRAN (R 3.6.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.6.0)
## mvabund * 4.1.3 2020-02-27 [1] CRAN (R 3.6.0)
## mvtnorm 1.0-11 2019-06-19 [1] CRAN (R 3.6.0)
## nlme * 3.1-141 2019-08-01 [1] CRAN (R 3.6.0)
## nnet 7.3-12 2016-02-02 [1] CRAN (R 3.6.1)
## openxlsx 4.1.4 2019-12-06 [1] CRAN (R 3.6.0)
## permute * 0.9-5 2019-03-12 [1] CRAN (R 3.6.0)
## pillar 1.4.3 2019-12-20 [1] CRAN (R 3.6.0)
## pkgbuild 1.0.6 2019-10-09 [1] CRAN (R 3.6.0)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.0)
## pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.6.0)
## plyr 1.8.6 2020-03-03 [1] CRAN (R 3.6.0)
## png 0.1-7 2013-12-03 [1] CRAN (R 3.6.0)
## pracma * 2.2.5 2019-04-09 [1] CRAN (R 3.6.0)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 3.6.0)
## processx 3.4.1 2019-07-18 [1] CRAN (R 3.6.0)
## ps 1.3.0 2018-12-21 [1] CRAN (R 3.6.0)
## pscl 1.5.5 2020-03-07 [1] CRAN (R 3.6.0)
## purrr * 0.3.4 2020-04-17 [1] CRAN (R 3.6.2)
## quantreg * 5.52 2019-11-09 [1] CRAN (R 3.6.0)
## R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.0)
## RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 3.6.0)
## Rcpp 1.0.4.6 2020-04-09 [1] CRAN (R 3.6.1)
## readr * 1.3.1 2018-12-21 [1] CRAN (R 3.6.0)
## readxl 1.3.1 2019-03-13 [1] CRAN (R 3.6.0)
## remotes 2.1.0 2019-06-24 [1] CRAN (R 3.6.0)
## reprex 0.3.0 2019-05-16 [1] CRAN (R 3.6.0)
## reshape 0.8.8 2018-10-23 [1] CRAN (R 3.6.0)
## rio 0.5.16 2018-11-26 [1] CRAN (R 3.6.0)
## rlang 0.4.6 2020-05-02 [1] CRAN (R 3.6.2)
## rmarkdown 2.1 2020-01-20 [1] CRAN (R 3.6.0)
## rpart 4.1-15 2019-04-12 [1] CRAN (R 3.6.1)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.0)
## rstudioapi 0.10 2019-03-19 [1] CRAN (R 3.6.0)
## rvcheck 0.1.8 2020-03-01 [1] CRAN (R 3.6.0)
## rvest 0.3.5 2019-11-08 [1] CRAN (R 3.6.0)
## scales 1.1.0 2019-11-18 [1] CRAN (R 3.6.0)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0)
## sp 1.3-2 2019-11-07 [1] CRAN (R 3.6.0)
## SparseM * 1.77 2017-04-23 [1] CRAN (R 3.6.0)
## statmod 1.4.34 2020-02-17 [1] CRAN (R 3.6.0)
## stringi 1.4.6 2020-02-17 [1] CRAN (R 3.6.0)
## stringr * 1.4.0 2019-02-10 [1] CRAN (R 3.6.0)
## survival 3.1-12 2020-04-10 [1] CRAN (R 3.6.2)
## testthat 2.2.1 2019-07-25 [1] CRAN (R 3.6.0)
## tibble * 3.0.1 2020-04-20 [1] CRAN (R 3.6.2)
## tidyr * 1.1.0 2020-05-20 [1] CRAN (R 3.6.2)
## tidyselect 1.1.0 2020-05-11 [1] CRAN (R 3.6.2)
## tidyverse * 1.3.0 2019-11-21 [1] CRAN (R 3.6.0)
## tweedie 2.3.2 2017-12-14 [1] CRAN (R 3.6.0)
## usethis 1.5.1 2019-07-04 [1] CRAN (R 3.6.0)
## vctrs 0.3.1 2020-06-05 [1] CRAN (R 3.6.2)
## vegan * 2.5-6 2019-09-01 [1] CRAN (R 3.6.0)
## weights 1.0.1 2020-02-12 [1] CRAN (R 3.6.0)
## withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.0)
## xfun 0.11 2019-11-12 [1] CRAN (R 3.6.0)
## xml2 1.3.0 2020-04-01 [1] CRAN (R 3.6.2)
## yaml 2.2.1 2020-02-01 [1] CRAN (R 3.6.0)
## zip 2.0.4 2019-09-01 [1] CRAN (R 3.6.0)
##
## [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library